﻿///@file Matrix.cpp
///@brief Implementation of the class representing a matrix
///@author Arnaud Duval
///@version 0.0
///@date 2010/02/26


#include <cmath>
#include <ctime>
#include <cstdlib>

#include "Matrix.h"

///@brief Default constructor
Matrix::Matrix()
{
	nbRows = 0;
	nbColumns = 0;
}

///@brief Constructor
///@param nrow number of rows
///@param ncol number of columns
Matrix::Matrix(unsigned int nrow, unsigned int ncol, double val)
{
	nbRows = nrow;
	nbColumns = ncol;
	arows = new	Vector[nbRows];
	for(unsigned int i = 0 ; i < nbRows ; i++)
	{
		arows[i] = Vector(nbColumns, val);
	}
}

///@brief Copy constructor
///@param m Matrix object to copy
Matrix::Matrix(const Matrix &m)
{
	nbRows = m.nbRows;
	nbColumns = m.nbColumns;
	arows = new Vector [nbRows];
	for(unsigned int i = 0 ; i < nbRows ; i++)
	{
		arows[i] = Vector (nbColumns);
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			arows[i][j] = m.arows[i][j];
	}
}

///@brief Copy constructor with temporary object
///@param m Temporary oject Matrix to copy
Matrix::Matrix(Matrix* m)
{
	nbRows = (*m).nbRows;
	nbColumns = (*m).nbColumns;
	arows = new Vector [nbRows];
	for(unsigned int i = 0 ; i < nbRows ; i++)
	{
		arows[i] = Vector (nbColumns);
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			arows[i][j] = (*m).arows[i][j];
	}
}

///@brief constructor of an special square matrix
///@param n size of square matrix
///@param c type of special matrix : I for identity matrix,
///R for random matrix,
///r for random symmetric
///D for operator that allow to compute the deviator of a stress Vector (works only for 3x3 matrix)
///P for operator that convert a stress vector to a strain vector (works only for 3x3 matrix)
Matrix::Matrix(unsigned int n, char c)
{
	switch(c)
	{
		case 'I':
		{
			nbRows = n;
			nbColumns = n;
			arows = new Vector [nbRows];
			for(unsigned int i = 0 ; i < nbRows ; i++)
			{
				arows[i] = Vector (nbColumns);
				for(unsigned int j = 1 ; j <= nbColumns ; j++)
				{
					if(i+1 ==j)	
					{
						arows[i][j] = 1.;
					}
					else
					{
						arows[i][j] = 0.;
					}
				}
			}
		}break;
		case 'R':
		{
			nbRows = n;
			nbColumns = n;
			arows = new Vector [nbRows];
			for(unsigned int i = 0 ; i < nbRows ; i++)
			{
				arows[i] = Vector (nbColumns);
				for(unsigned int j = 1 ; j <= nbColumns ; j++)
				{
					arows[i][j] = rand()/(double)RAND_MAX;
				}
			}
		}break;
		case 'r':
		{
			nbRows = n;
			nbColumns = n;
			arows = new Vector[nbRows];
			for(unsigned int i = 0 ; i < nbRows ; i++)
				arows[i] = Vector(nbColumns);
			for(unsigned int i = 0 ; i < nbRows ; i++)
			{
				for(unsigned int j = i+1 ; j <= nbColumns ; j++)
				{
					arows[i][j] = rand()/(double)RAND_MAX;
					if(j != i+1)
					{
						arows[j-1][i+1] = arows[i][j];
					}
				}
			}
		}break;
        case 'D':
        {
            assert(n == 3);
            nbRows = n;
            nbColumns = n;
            arows = new Vector[nbRows];
            for(unsigned int i = 0 ; i < nbRows ; i++)
                arows[i] = Vector(nbColumns);
            arows[0][1] = 2./3.;
            arows[0][2] = -1./3.;
            arows[0][3] = 0.;
            arows[1][1] = -1./3.;
            arows[1][2] = 2./3.;
            arows[1][3] = 0.;
            arows[2][1] = 0.;
            arows[2][2] = 0.;
            arows[2][3] = 1.;
        }break;
        case 'P':
        {
            assert(n == 3);
            nbRows = n;
            nbColumns = n;
            arows = new Vector[nbRows];
            for(unsigned int i = 0 ; i < nbRows ; i++)
                arows[i] = Vector(nbColumns);
            arows[0][1] = 1.;
            arows[0][2] = 0.;
            arows[0][3] = 0.;
            arows[1][1] = 0.;
            arows[1][2] = 1.;
            arows[1][3] = 0.;
            arows[2][1] = 0.;
            arows[2][2] = 0.;
            arows[2][3] = 2.;
        }break;
		default:
		{}
	}
}

///@brief Destructor
Matrix::~Matrix()
{
	if( nbRows != 0)
		delete[] arows;
}

///@brief Overload of () operator, direct access to a component
///@param r index of raw
///@param c index of column
///@return Direct access to a component
double& Matrix::operator()(unsigned int r, unsigned int c)
{
	assert((c > 0)&&(r > 0));
	assert ((r <= nbRows)&&(c <= nbColumns));
	return arows[r - 1][c];
}

///@brief Overload of [] operator : return a row
///@param i index of row to return
///@return a row
Vector& Matrix::operator[](unsigned int i) const
{
	assert((i > 0)&&(i <= nbRows));
	return arows[i - 1];
}	


///@brief Overload of = operator
///@param m Right hand operand
Matrix& Matrix::operator = (const Matrix& m)
{
	if((nbColumns == 0)||(nbRows == 0))
	{
		nbRows = m.nbRows;
		nbColumns = m.nbColumns;
		arows = new Vector [nbRows];
		for(unsigned int i = 0 ; i < nbRows ; i++)
			arows[i] = Vector(nbColumns);
	}
	else
	{
		assert((nbColumns == m.nbColumns)&&(nbRows == m.nbRows));
	}

	for(unsigned int i = 0 ; i < nbRows ; i++)
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			arows[i][j] = m.arows[i][j];

	return (*this);
}
		
///@brief Overload of = operator for temporary objects
///@param m Right hand operator
Matrix& Matrix::operator = (Matrix* m)
{
	if((nbColumns == 0)||(nbRows == 0))
	{
		nbRows = (*m).nbRows;
		nbColumns = (*m).nbColumns;
		arows = new Vector [nbRows];
		for(unsigned int i = 0 ; i < nbRows ; i++)
			arows[i] = Vector(nbColumns);
	}
	else
	{
		assert((nbColumns == (*m).nbColumns)&&(nbRows == (*m).nbColumns));
	}

	for(unsigned int i = 0 ; i < nbRows ; i++)
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			arows[i][j] = (*m).arows[i][j];

	return (*this);
}

///@brief Copies Matrix without size check (use with caution)
///@param m Matrix object to copy
Matrix& Matrix::operator &= (const Matrix& m)
{
	if(!m)
	{
		if(nbRows)
			delete[] arows;
	}
	else
	{
		if(nbRows)
			delete[] arows;
		arows = new Vector[nbRows = m.nbRows];
		nbColumns = m.nbColumns;
		for(unsigned int i = 1 ; i <= m.nbRows; i++)
			(*this)[i] &= m[i];
	}
	return *this;
}

///@brief Overload of += operator
///@param m right hand operator
void Matrix::operator += (const Matrix& m)
{
	assert(nbColumns == m.nbColumns);
	assert(nbRows == m.nbRows);
	for(unsigned int i = 0 ; i < nbRows ; i++)
		arows[i] += m.arows[i];
}
///@brief Overload of -= operator
///@param m right hand operator
void Matrix::operator -= (const Matrix& m)
{
	assert(nbColumns == m.nbColumns);
	assert(nbRows == m.nbRows);
	for(unsigned int i = 0 ; i < nbRows ; i++)
		arows[i] -= m.arows[i];
}

///@brief return number of rows
unsigned int Matrix::NbRows() const
{
	return nbRows;
}

///@brief return nomber of columns
unsigned int Matrix::NbColumns() const
{
	return nbColumns;
}


///@brief extract a row from the matrix
///@param i index of the row to extract
Vector Matrix::ExtractRow(unsigned int i)
{
	assert((i > 0)&&(i<=nbRows));
	return arows[i-1];
}

///@brief extract a column from the matrix
///@param j index of the column to extract
Vector Matrix::ExtractColumn(unsigned int j)
{
	assert((j > 0)&&(j <= nbColumns));
	Vector res(nbRows);
	for(unsigned int i = 1 ; i <= nbRows ; i++)
		res[i] = (arows[i-1])[j];
	return res;
}

///@brief Extract diagonal from the matrix
Vector Matrix::ExtractDiag()
{
	assert(nbColumns == nbRows);
	Vector res(nbRows);
	for(unsigned int i = 1 ; i <= nbRows ; i++)
		res[i] = arows[i-1][i];
	return res;
}


///@brief swap two rows
///@param i1 index of first row
///@param i2 index of second row
void Matrix::swapRow(unsigned int i1, unsigned int i2)
{
	assert((i1 > 0)&&(i1<=nbRows));
	assert((i2 > 0)&&(i2<=nbRows));
	Vector temp;
	temp = arows[i1-1];
	arows[i1-1] = arows[i2-1];
	arows[i2-1] = temp;
}

///@brief swap two columns
///@param i1 index of first column
///@param i2 index of second columns
void Matrix::swapColumn(unsigned int i1, unsigned int i2)
{
	assert((i1 > 0)&&(i1<=nbColumns));
	assert((i2 > 0)&&(i2<=nbColumns));
	double temp;
	for(unsigned int i = 0 ; i < nbRows ; i++)
	{
		temp = (arows[i])[i1];
		(arows[i])[i1] = (arows[i])[i2];
		(arows[i])[i2] = temp;
	}
}

///@brief Invert matrix
Matrix Matrix::Invert()
{
	return InvertGaussJordan();
}

///@brief Invert Matrix with Gauss-Jordan method
Matrix Matrix::InvertGaussJordan()
{
	///@brief Gauss-Jordan method
	assert(nbRows == nbColumns);
	Matrix temp = *this;
	Matrix res(nbRows, 'I');
	
	double a,b;
	a = 0.;
	b = 0.;
	int c = 0;
	
	for(unsigned int k = 1 ; k <= nbRows ; k++)
	{
		///@brief Pivot determination
		a = temp(k,k);
		c = 0;
		while((std::abs(a) < 1.E-30)&&(k+c < this->NbRows()))
		{
			c++;
			for(unsigned int q = 1 ; q <= nbRows ; q++)
			{
				temp(k,q) = temp(k,q) + temp(k+c, q);
				res(k,q) = res(k,q) + res(k+c,q);
			}
			a = temp(k,k);
		}
		
		///@brief rows normalization
		
		for(unsigned int l = 1 ; l <= nbRows ; l++)
		{
			temp(k,l) = temp(k,l)/a;
			res(k,l) = res(k,l)/a;			
		}
		
		
		///@brief Gauss Jordan reduction
		
		for(unsigned int i = 1; i <= nbRows ; i++)
		{
			b = temp(i,k);
			if(i != k)
			{
				for(unsigned int j = 1 ; j <= nbRows ; j++)
				{
					temp(i,j) = temp(i,j) - b*temp(k,j);
					res(i,j) = res(i,j) - b*res(k,j);
				}
			}
		}
	}
	return res;
}




///@brief Invert matrix with Cholesky method (for symmetric and positive-defined matrices)
Matrix Matrix::InvertCholesky()
{
	///@brief Cholesky
	assert(nbRows == nbColumns);
	double ztol = 1.e-10;
	Matrix t(nbRows, nbRows, 0.);

    for(unsigned int i = 0 ; i < nbRows ; i++)
    {
        double s = 0.;
        for(unsigned int k = 0 ; k < i ; k++)
            s = s + t.arows[k][i+1]*t.arows[k][i+1];
        
        double d = arows[i][i+1] - s;
        if(fabs(d) < ztol)
            t.arows[i][i+1] = 0.;
        else
            if (d < 0.)
                std::cerr << "Matrix not positive-defined";
            t.arows[i][i+1] = sqrt(d);
    
        for(unsigned int j = i+1 ; j < nbRows ; j++)
        {
            s = 0.;
            for(unsigned int k = 0 ; k < i ; k++)
                s = s + t.arows[k][i+1] * t.arows[k][j+1];
                
            if(fabs(s) < ztol)
                s = 0.0;
            t.arows[i][j+1] = (arows[i][j+1] - s)/t.arows[i][i+1];
        }
    }

	Matrix res(nbRows, nbRows, 0.);

    // Backward step for inversion
    for(int j = nbRows-1 ; j >= 0 ; j--)
    {
        double tjj = t.arows[j][j+1];
        double s = 0.;
        for(unsigned int k = j+1 ; k < nbRows ; k++)
            s = s + t.arows[j][k+1]*res.arows[j][k+1];
        res.arows[j][j+1] = 1./(tjj*tjj) - s/tjj;
        for(int i = j-1 ; i >= 0 ; i--)
        {
            double sum = 0.;
            for(unsigned int k = i+1 ; k < nbRows ; k++)
            {
                sum = sum + t.arows[i][k+1]*t.arows[k][j+1];
            }
            res.arows[j][i+1] = res.arows[i][j+1] = -sum/t.arows[i][i+1];
        }
    }


	return res;
}


///@brief Invert matrix with analytic method (limited to 2x2 and 3x3 matrices)
Matrix Matrix::InvertAnalytic()
{
	assert(nbColumns == nbRows);
	assert((nbRows == 2)||(nbRows == 3));

	Matrix res(nbRows, nbRows);

	if(nbRows == 2)
	{
		double det = arows[0][1]*arows[1][2] - arows[0][2]*arows[1][1];
		res.arows[0][1] = arows[1][2]/det;
		res.arows[0][2] = -arows[0][2]/det;
		res.arows[1][1] = -arows[1][1]/det;
		res.arows[1][2] = arows[0][1]/det;
	}
	else
	{
		double det = arows[0][1]*(arows[1][2]*arows[2][3] - arows[1][3]*arows[2][2])
                    -arows[0][2]*(arows[1][1]*arows[2][3] - arows[1][3]*arows[2][1])
					+arows[0][3]*(arows[1][1]*arows[2][2] - arows[1][2]*arows[2][1]);
		res.arows[0][1] = (arows[1][2]*arows[2][3] - arows[1][3]*arows[2][2])/det;
		res.arows[0][2] = -(arows[0][2]*arows[2][3] - arows[0][3]*arows[2][2])/det;
		res.arows[0][3] = (arows[0][2]*arows[1][3] - arows[0][3]*arows[1][2])/det;
		res.arows[1][1] = -(arows[1][1]*arows[2][3] - arows[1][3]*arows[2][1])/det;
		res.arows[1][2] = (arows[0][1]*arows[2][3] - arows[0][3]*arows[2][1])/det;
		res.arows[1][3] = -(arows[0][1]*arows[1][3] - arows[0][3]*arows[1][1])/det;
		res.arows[2][1] = (arows[1][1]*arows[2][2] - arows[1][2]*arows[2][1])/det;
		res.arows[2][2] = -(arows[0][1]*arows[2][2] - arows[0][2]*arows[2][1])/det;
		res.arows[2][3] = (arows[0][1]*arows[1][2] - arows[0][2]*arows[1][1])/det;
	}

	return res;
}

///@brief Transpose matrix
Matrix Matrix::Transpose()
{
	assert((nbColumns != 0)&&(nbRows != 0));
	Matrix res(nbColumns, nbRows);
	for(unsigned int i = 0 ; i < nbRows ; i++)
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			res(j, i+1) = (arows[i])[j];
	return res;
}

///@brief Compute matrix determinant
double Matrix::Det()
{
	assert(nbColumns != 0);
	///@brief Matrix must be squared
	assert(nbColumns == nbRows);
	double res = 0.;
	
	if (nbRows==1) 
		return (arows[0])[1];
	
	Matrix m(nbRows-1, nbColumns-1);

	for(unsigned int i = 1 ; i <= nbRows ; i++)
	{
		for(unsigned int j = 2 ; j <= nbColumns ; j++)
			for(unsigned int k=1, l=1 ; k <= nbRows ; k++)
				if(k != i)
					m(l++, j-1)	 = (arows[k-1])[j];
		
		if((arows[i-1])[1])
		{
			if( ((float)(i)/2) == int((float) (i)/2) )
				res -= (arows[i-1])[1] * m.Det();
			else
				res += (arows[i-1])[1] * m.Det();
		}
		
	}
	return res;
}

///@brief Compute Frobenius norm
double Matrix::FrobeniusNorm()
{
	assert(nbColumns != 0);
	///@brief Matrix must be squared
	assert(nbColumns == nbRows);
	double res = 0.;
	for(unsigned int i = 1 ; i <= nbColumns ; i++)
		for(unsigned int j = 1 ; j <= nbColumns ; j++)
			res += pow((arows[i-1])[j], 2.);

	res = sqrt(res);
	return res;
}

///@brief Add a submatrix into current matrix
///@param mat Submatrix to be inserted
///@param r Index of the row where 1,1 component of submatrix will be inserted
///@param c Index of the column where 1,1 component of submatrix will be inserted
void Matrix::AddSubMatrix(Matrix mat, unsigned int r, unsigned int c)
{
	for(unsigned int i = 1 ; i <= mat.nbRows ; i++)
		for(unsigned int j = 1 ; j <= mat.nbColumns ; j++)
		{
			arows[i+r-2][j+c-1] += mat(i,j);
		}
}

///@brief Compute QR decomposition
///@param Q Orthogonal Q result Matrix
///@param R triangular sup R result Matrix
///@warning Must be tested for rectangular matrices
void Matrix::QRDecomposition(Matrix &Q, Matrix &R)
{
	unsigned int n = nbRows;
	unsigned int m = nbColumns;

	Matrix H(m,'I');
	R = this;
	Vector v(m);
	for(unsigned int k = 1 ; k <= n-1 ; k++)
	{
		double alpha = 0.;
		for(unsigned int i = k ; i <= m ; i++)
			alpha += pow(R(i,k), 2.);
		alpha = sqrt(alpha);
		double beta = pow(alpha, 2.) - alpha*R(k,k);
		v[k] = R(k,k)-alpha;
		for(unsigned int i = k+1 ; i <= m ; i++)
			v[i] = R(i,k);
		for(unsigned j = k ; j <= n ; j++)
		{
			double c = 0.;
			for(unsigned int i = k ; i <= m ; i++)
				c+= v[i]*R(i,j);
			c *= 1./beta;
			for(unsigned int i = k ; i <= m ; i++)
				R(i,j) = R(i,j) - c*v[i];
		}
		for(unsigned int j = 1 ; j <= m ; j++)
		{
			double c = 0.;
			for(unsigned int i = k ; i <= m ; i++)
				c+= v[i]*H(i,j);
			c *= 1./beta;
			for(unsigned int i = k ; i <= m ;i++)
				H(i,j) = H(i,j) - c* v[i];
		}
	}
	Q = H.Transpose();
}

///@brief Compute Lanczos algorithm for eigenvalues extraction
///@param T Tridiagonal output matrix
///@param V Base change matrix
void Matrix::Lanczos(Matrix &T, Matrix &V)
{
	assert(nbRows == nbColumns);
	Vector *v = new Vector[nbColumns+2];
	Vector *w = new Vector[nbColumns+2];
	Matrix A = this;
	v[0] = Vector(nbColumns);
	v[1] = Vector(nbColumns);
	Vector beta(nbColumns+1);
	Vector alpha(nbColumns+1);
	(v[1])[1] = 1.;
	beta[1] = 0.;
	unsigned int m = nbColumns;
	for(unsigned int j = 1 ; j <= m ; j++)
	{
		w[j] = (A|v[j]) - beta[j]*v[j-1];
		alpha[j] = w[j]|v[j];
		w[j] = w[j] - alpha[j]*v[j];
		beta[j+1] = (w[j]).Euclide_Norm();
		v[j+1] = w[j]/beta[j+1];
	}
	T = Matrix(m,m);
	V = Matrix(m,m);
	for(unsigned int i = 1 ; i <= m ; i++)
	{
		T(i,i) = alpha[i];
		if(i != nbColumns)
		{
			T(i,i+1) = beta[i+1];
			T(i+1,i) = beta[i+1];
		}
		for(unsigned int j = 1 ; j <= m ; j++)
		{
			V(i,j) = (v[j])[i];
		}
	}
}

///@brief Check if Matrix is symmetric, assuming a tolerance
///@param toler tolerance to decide if two components are equal
bool Matrix::IsSymmetric(double toler)
{
	if(nbColumns != nbRows)
		return false;
	for(unsigned int i = 1 ; i <= nbColumns ; i++)
	{
		for(unsigned int j = i+1 ; j <= nbRows ; j++)
		{
			if (fabs(arows[i-1][j] - arows[j-1][i]) > toler)
				return false;
		}
	}
	return true;
}

///@brief Set all components to zero
void Matrix::Nullify()
{
  
    for(unsigned int i = 1 ; i <= nbColumns ; i++)
        for(unsigned int j = 1 ; j <= nbRows ; j++)
            arows[j-1][i] = 0.;
}

///@brief Return true if Matrix exists
///@param m Matrix object to be tested
bool operator ! (const Matrix& m)
{
	return !m.nbRows;
}

///@brief Overload of + operator
///@param a left hand operand
///@param b right hand operand
Matrix operator+ (const Matrix& a, const Matrix& b)
{
	assert(a.nbRows == b.nbRows);
	assert(a.nbColumns == b.nbColumns);
	Matrix res(a.nbRows, a.nbColumns);
	for(unsigned int i = 1 ; i <= a.nbRows ; i++)
		for(unsigned int j = 1 ; j <= a.nbColumns ; j++)
			res.arows[i-1][j] = a.arows[i-1][j] + b.arows[i-1][j];
	return res;
}
		
		
///@brief Overload of - operator
///@param a left hand operand
///@param b right hand operand
Matrix operator- (const Matrix& a, const Matrix& b)
{
	assert(a.nbRows == b.nbRows);
	assert(a.nbColumns == b.nbColumns);
	Matrix res(a.nbRows, a.nbColumns);
	for(unsigned int i = 1 ; i <= a.nbRows ; i++)
		for(unsigned int j = 1 ; j <= a.nbColumns ; j++)
			res.arows[i-1][j] = a.arows[i-1][j] - b.arows[i-1][j];
	return res;
}

///@brief Overload of * operator, cross product of two vector A_ij = B_i * C_j
///@param v1 first vector
///@param v2 second vector
Matrix operator * (const Vector &v1, const Vector &v2)
{
	assert((v1.Dim() > 0)&&(v2.Dim() > 0));
	Matrix res(v1.Dim(), v2.Dim());
	for(unsigned int i = 1 ; i <= v1.Dim() ; i++)
		for(unsigned int j = 1 ; j <= v2.Dim() ; j++)
		{
			res.arows[i-1][j]	= v1[i]*v2[j];
		}
	return res;
}
		
		
///@brief overload of | operator, contracted product : A_i = B_ij * C_j
///@param m matrix
///@param v vector
Vector operator| (const Matrix& m, const Vector& v)
{
	assert(v.Dim() > 0);
	assert(m.NbColumns() == v.Dim());
	assert(m.NbRows() > 0);
	Vector res(m.NbRows());
	for(unsigned int i = 1 ; i <= m.NbRows() ; i++)
	{
		res[i] = 0.;
		for(unsigned int j = 1 ; j <= m.NbColumns() ; j++)
		{
			res[i] += m.arows[i-1][j] * v[j];
		}
	}
	return res;
}

///@brief overload of | operator, contracted product : A_i = c_j * B_ij
///@param v vector
///@param m matrix
Vector operator| (const Vector& v, const Matrix& m)
{
	assert(v.Dim() > 0);
	assert(m.NbColumns() == v.Dim());
	assert(m.NbRows() > 0);
	Vector res(m.NbRows());
	for(unsigned int i = 1 ; i <= m.NbRows() ; i++)
	{
		res[i] = 0.;
		for(unsigned int j = 1 ; j <= m.NbColumns() ; j++)
		{
			res[i] += m.arows[i-1][j] * v[j];
		}
	}
	return res;
}

///@brief Overload of || operator : twice contracted product A = B_ij * C_ij
///@param B matrix B_ij
///@param C matrix C_ij
double operator|| (const Matrix& B, const Matrix& C)
{
	assert(B.nbColumns == C.nbColumns);
	assert(B.nbRows == C.nbRows);
	assert(B.nbRows == B.nbColumns);
	assert(B.nbColumns > 0);
	double A = 0.;
	for(unsigned int i = 1 ; i <= B.nbRows ; i++)
	{
		for(unsigned int j = 1 ; j <= B.nbColumns ; j++)
		{
			A += B.arows[i-1][j] * C.arows[i-1][j];
		}
	}
	return A;
}

///@brief Overload of | operator : contracter product : A_ij = B_ik * C_kj
///@param a left hand operand
///@param b right hand operand
Matrix operator | (const Matrix& a, const Matrix& b)
{
	assert(a.nbColumns == b.nbRows);
	
	Matrix res(a.nbRows, b.nbColumns);
	
	for(unsigned int i = 1; i <= a.nbRows ; i++)
		for(unsigned int j = 1 ; j <= b.nbColumns ; j++)
		{
			res(i,j) = 0.;
			for(unsigned int k = 1 ; k <= a.nbColumns ; k++)
				res(i,j) += a.arows[i-1][k] * b.arows[k -1][j];
		}
	
	return res;
}

///@brief Overload of * operator : product scalar * matrix
///@param s scalar
///@param m matrix
Matrix operator * (const double& s, const Matrix &m)
{
	Matrix res(m.nbRows, m.nbColumns);
	for(unsigned int i = 0 ; i < m.nbRows ; i++)
		res.arows[i] = s*m.arows[i];
	return res;
}
		
///@brief Overload of * operator : product matrix * scalar
///@param m matrix
///@param s scalar
Matrix operator * (const Matrix& m, const double &s)
{
	Matrix res(m.nbRows, m.nbColumns);
	for(unsigned int i = 0 ; i < m.nbRows ; i++)
		res.arows[i] = s*m.arows[i];
	return res;
}

///@brief Overload of / operator : division Matrix / scalar
///@param m left hand operand : Matrix
///@param s right hand operand : scalar
Matrix operator / (const Matrix& m, const double& s)
{
	assert(s != 0.);
	Matrix res(m.nbRows, m.nbColumns);
	for(unsigned int i = 0 ; i < m.nbRows ; i++)
		res.arows[i] = m.arows[i] / s;
	return res;
}

///@brief overload of << operator
///@param f out stream
///@param m Matrix to put in the out stream
std::ostream& operator<<(std::ostream &f, const Matrix& m)
{
	for(unsigned int i = 0 ; i < m.nbRows ; i++)
	{
		f << m.arows[i];
		if(i != (m.nbRows - 1))
		{
			f << "\n";
		}
	}
	
	return f;
}